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The spontaneous symmetry breaking taking place in the direction perpendicular to the energy 
fiux in a dilute vibrofluidized granular system is investigated, using both a hydrodynamic description 
and simulation methods. The latter include molecular dynamics and direct Monte Carlo simulation 
of the Boltzmann equation. A marginal stability analysis of the hydrodynamic equations, carried 
out in the WKB approximation, is shown to be in good agreement with the simulation results. The 
shape of the hydrodynamic profiles beyond the bifurcation is discussed. 
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I. INTRODUCTION 

Granular materials are assemblies of macroscopic particles dissipating their energy through inelastic coUisions [0. 
They exhibit a very rich phenomenology that is only partially understood. One of the most peculiar behaviors of 
granular systems, which has attracted a lot of attention in the last years, is their tendency to spontaneously develop 
strong spatial inhomogeneities. In many different situations, the density shows a sharp profile that is not induced by 
the boundary conditions. This phenomenon is often referred to as a clustering effect since high density regions 
coexist in the system with regions where the density is very low. 

In vibrated granular systems, clustering effects show up in many cases as a spontaneous symmetry breaking in the 
direction parallel to the vibrating wall. Consider a gas enclosed in a box that is being supplied energy through a 
vibrating wall located at x = 0. There are no other external forces acting on the system. The box is divided into 
two equal compartments by a wall normal to the x axis starting at a certain height. At sufhciently low average 
density, the hydrodynamic fields are symmetric to both sides of the partition, but above a critical average density, 
which depends on the value of the restitution coefficient, an asymmetry in the number of particles at each side of 
the container occurs |^ . This asymmetry has been shown to be associated with a bifurcation of the solution of the 
hydrodynamic equations describing the state of the system. 

A similar symmetry breaking has been observed in a system in presence of a gravitational force acting in the x 
direction. In this case, the system is unbounded for a; > 0, and the partition has a hole at a certain height. Again, 
an asymmetry in the number of particles in the two compartments develops if a control parameter, dependent on the 
amplitude of the vibration and the degree of inelasticity, is larger than a critical value ^ ,|^ . 

Symmetry breaking in the direction parallel to the vibrating wall has also been observed in systems without any 
partition of the container. Sunthar and Kumaran |^ have reported molecular dynamics simulation results showing 
the presence of convection rolls and phase separation into coexisting dense and dilute regions in a granular system in 
presence of gravity. The phase separation takes place on the surface of the vibrating wall and, as already indicated, 
in the direction perpendicular to the energy flux. No theoretical explanation for this phenomenon is provided in ref. 

, although the effect of the different parameters controlling the behavior of the system is discussed in detail, on the 
basis of the simulation results. For a two-dimensional closed system in absence of gravity, a transversal continuous 
spontaneous symmetry breaking has also been predicted |^]. Not at all surprisingly, the gradients are now sharper 
next to the elastic wall, opposite to the energy source. This is consistent with the positions that must have the 
holes in systems with a separating wall in order to observe the symmetry breaking, with and without gravitational 
field acting on the system. The work by Livne et al. is restricted to the nearly elastic limit and it is based on a 
numerical marginal stability analysis of the hydrodynamic equations. The predictions of this analysis are compared 
with numerical solutions of the own hydrodynamic equations with the appropriate boundary conditions. 

In this paper, the bifurcation predicted in ref. will be considered again. There are several reasons for that. 
Firstly, hydrodynamic equations derived from the Boltzmann equation for smooth inelastic hard disks and valid, in 
principle, for arbitrary inelasticity will be used, then somehow extending the previous results. Secondly, instead of a 
numerical analysis of the stability of the solutions of the hydrodynamic equations, an analytical study, based on the 
WKB approximation, will be presented here. One of the main advantages of this approach is that the dimensionless 
control parameter governing the bifurcation phenomenon is clearly identified. A third motivation for the present work 
is to report simulation results, both from molecular dynamics and also from Monte Carlo simulation of the Boltzmann 
equation, showing the existence of the predicted transition. Since these simulation techniques do not contain any 
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externally introduced hydrodynamic concept, they provide a direct proof of the existence of the continuous symmetry 
breaking, and a test of the theoretical predictions. Attention will be also paid to the form of the hydrodynamic profiles 
beyond the bifurcation. This leads to a deeper understanding of the development of the instability. In any case, it is 
clear that the work in ref. opened the way to more systematic investigations of this instability, like the one in this 
paper. 

The plan of the paper is as follows. In Sec. ||, the hydrodynamic description of the one-dimensional state of a low 
density vibrofluidized granular gas will be briefly summarized. The analytical expressions of the hydrodynamic profiles 



are given. This state is the starting point for the marginal stability analysis developed in Sec. [11. Linearization of the 
hydrodynamic equations around the one-dimensional state leads to a second order linear differential equation. The 
WKB solution of the closed problem posed by this equation and the corresponding boundary conditions is built up. 
This requires to consider three different cases, depending on the value of the parameters characterizing the system. 
From the WKB solution, the marginal stability curve follows easily. Simulation results are presented and compared 
with the theoretical predictions in Sec. where an order parameter characterizing the transition is defined. A good 
agreement is found. The last section contains some final remarks, as well as the comparison of the results derived 
here with those in ref. in the common range of applicability, namely nearly elastic collisions and very low density. 

II. BASIC EQUATIONS AND THE REFERENCE STATE 

For a steady state without macroscopic flows, the balance equations for a two-dimensional gas of smooth inelastic 
hard disks of mass m and diameter a have the form 

V • P = 0, (1) 

V-q + nTC = 0, (2) 

where n and T are the number density and granular temperature (with Boltzmann's constant set equal to unity), 
respectively. In the low density limit, for a gas described by the inelastic Boltzmann equation, and to lowest order in 
the gradients (Navier-Stokes order), the pressure tensor P and heat flux q, for the steady state under consideration 
are given by ||,§ 

P=pl, (3) 

q = -kVT - /iVn, (4) 

I being the unit tensor, p = nT the hydrodynamic pressure, k the heat conductivity, and /i a transport coefficient 
that has no analogue in the elastic case. These transport coefficients are proportional to the elastic heat conductivity 

K ^ K*ia)Ko{T), fi = y*{a.) ^'^°^^\ (5) 



2 / T 



1/2 

^o{T) = -(—)■ (6) 
a \TTm J 

Finally, the cooling rate in the same approximation is related with the elastic shear viscosity r/o by 

c = r(«)-, (7) 

Vo 

, , 1 /mr\'/' 

^°(^)-2^(-) ■ 

The functions k*, /i*, and C,* only depend on the constant coefficient of normal restitution a characterizing the 
inelasticity of collisions. Their explicit expressions are given in refs. and |10|, and will be not reproduced here. By 
using the above expressions we get from Eqs. dll) and (0), 
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Next, we specify the boundary conditions. We consider that the system has TV particles enclosed in a rectangular box 
with dimensions and Ly. The wall located at x = is vibrating and, therefore, supplies energy to the system. 
This energy is needed in order to keep and sustain a fluidized steady state. The other three walls, located at a; = La;, 
y = Q, and y = Ly, respectively, are at rest. For the sake of simplicity, collisions of particles with all four walls are 
assumed to be elastic. Then, the mathematical boundary conditions to be imposed are: 



dT 
dx 



dT 

dyjy=. 



dT 
dyjy^i 



(11) 



[K*{a) — ^J,*{a)] 



dT' 



= Q. 



(12) 



x=0 



Equations ( pT| ) express that the heat flux must vanish at the inmobile walls, while Eq. (jl^) is the energy balance at 
the vibrating wall. The quantity Q is the rate of energy input through this wall per unit of length. Its calculation in 
terms of the parameters defining the motion of the wall has been addressed in several works. Here we will consider 
the simplest possibility, namely that the wall moves in a sawtooth manner with velocity vj, [^-12|. This is a good 
approximation to more realistic motions, as long as the characteristic frequency of vibration of the wall be much larger 
than the collision rate of the gas molecules in its vicinity. In addition, it will be assumed that the amplitude of the 
vibration is much smaller than the mean free path of the particles of the gas next to it, so collective motions in the 
system are not being generated. Under the above conditions, it is 



Q = pvb- 



(13) 



The closed mathematical problem defined by Eqs. ([9|)-([l^) admits a y-independent solution T = Tb,{x) that has 
been discussed in detail in rcf. pO[ |. The existence of this one-dimensional solution had been previously noticed by 
Grossman et al. [|l^. In the following, we will refer to this solution as the reference state and its properties will be 
characterized with a subindex R. The reference temperature profile is 



Tr{£,x) — To,R 
where is a dimensionless scaled length defined by 



cosh(g* - £,x) 
cosh ^* 



with \r{x) being the local mean free path. 



and 



dx 



^r{x) 
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Moreover, ^* is the value of for x = L^, i.e.. 



2^/2anR{x) 
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(14) 



(15) 



(16) 



(17) 



(18) 



Therefore, ^* is proportional to the number of monolayers of particles perpendicular to the x-axis at rest, crN/Ly. 
Finally, the uniform pressure of the reference state is 
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To,fi(2e*+sinh2g) 
8y/2a{a)aLx cosh^ ^* 

and the temperature Tq./j of the gas next to the vibrating wall is given by 



(19) 



The simplicity of the above results is a consequence of the limiting kind of motion of the vibrating wall considered, 
whose only only effect is to transfer energy to the grains, without inducing any periodic motion in the system. In some 
previous studies [^Q, a "thermal" wall, instead of a vibrating one, has been considered at x 0. By definition, 
particles which collide with a thermal wall leave it with the velocity distribution corresponding to the temperature 
of the wall. Although this kind of walls are far from reality for granular systems, their consideration might be useful 
for comparison purposes. The only change to be made in the above discussion in order to apply it to a system with 
a thermal wall, is to replace the boundary condition given in Eq. (^) by the requirement that the temperature Tq.r 
has the value determined by the wall. As a consequence, the expressions of the hydrodynamic profiles remain the 
same, while Eq. ( po|) does not apply in this case. 

III. MARGINAL STABILITY ANALYSIS 

Our aim now is to investigate whether the system described in the previous Section exhibits another steady state, 
in addition to the reference one. Then, we introduce a perturbation dT(x, y) by 

T{x,v)^TR(x)+5T{x,y), (21) 

with ST{x, y) <C Tii{x), and search for a solution of this form to Eqs. ([9|)-(p^ with the boundary conditions ( pT| ) and 
( |l^ ) . Substitution of Eq. (|2l| ) into Eq. ( |lo|) and hnearization in ST yields 



92 ^^lfdlnTR\^ , fXR\^ 52 



ST = —Sp, (22) 



where Sp has been defined hy p ~ pii + Sp, is the dimensionless scale introduced in Eq. (llSD , and 



Here A is the average mean free path. 



= Vaiaj | . (23) 



A = , (24) 

n — N/LxLy. Next, we consider factorized solutions of Eq. (p^), 

ST{^x,iy)^HM^v)- (25) 
The conditions at the boundaries to be verified by the functions cj) and ip are 



d^J^^=i^^ V^ej^.^o y^^yJ^y^i; 



where 
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ey=2^/M^—. (28) 
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Moreover, Eqs. (|11|) and Jll 
perturbation. When Eq. (12! 



also imply that it must be 5p = 0, i.e. the pressure is not changed by the small 
is substituted into Eq. (^2|), use of separation of variables leads to the equations 



1 d^^iCy) 



da 



(29) 



1 f dlnTg 
dQ 2 di. 
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(30) 



where k is the constant of separation. The solution of Eq. (29) verifying the corresponding boundary conditions in 
Eq. <M) is 



ip — A cos fc^j,, 



(31) 



with A an arbitrary constant. Moreover, the values of k are restricted to fc = nq/^*, q being an integer. It is important 
to realize that the low density limit does not imply either ^* ^ 1 or ^* <C 1. In fact, both quantities can be large in 
a very dilute system. An additional restriction to be required to 5T{x, y) is that the total number of particles in the 
system, iV, is conserved. It is easily seen that this condition is equivalent to 



d^x I ^6 



(32) 



Equation ( pSj) with ip{£,y) given by Eq. (^) guarantees that this equality is verified. 

An equation having a structure similar to Eq. ( |30| ) was obtained in ref. for a dense system in the limit of 
nearly elastic collisions, by employing approximate constitutive relations introduced by Grossman et al. [ p^ . While 
the equation in [pj was solved using numerical techniques, here we will use a Wentzel, Kramers, Brillouin (WKB) 
approximation to investigate the possible solutions ||l^. First we use Eqs. ( p^ ) and (19) to rewrite Eq. ( po| ) as 



dil 



with 



le^c'cosh^CC^-C.) 



COSh'(^J -Cx) 



(2^2 +sinh 2^5)2 



(33) 



(34) 



which is a monotonic increasing function of in all the interval < l£ 

To construct the WKB exponential approximation of Eq. (|3^) it is necessary to consider three different ranges of 
parameters, that will be discussed separately in the following. 

a) The function f{£,x) is positive everywhere in the system. This is equivalent to /(O) > or 



k < 



2g + sinh2C; 
2V2eScosh3eS 



Then, the WKB solution is oscillatory, 

m) = 



Cl 



exp 



d£.'x9{Q 



C2 



V5(C^ 

where ci and C2 are constants, and 

Imposing the boundary condition at = Cx leads to 



: CXp 



d^'x9{0 



(35) 



(36) 



(37) 



(38) 
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with c another constant. When the boundary condition a,t = 0, Eq. (p7|), is also required, the consistency condition 

3tanh^* 



tan 



<73(0)cosh2ej 



(39) 



follows. The above equation determines the possible values of the parameters for which a WKB solution of the 
differential equation (|3^) exists in the region under study. 

If a thermal wall is considered at a; = 0, the boundary condition ( p7| ) must be replaced by (/)(0) — 0, and instead of 
Eq. it is found: 



I.e. 



cos / d^x9{D = 0^ 
/o 



(40) 



(41) 



where q is an arbitrary integer. 

b) The function f{^x) is negative everywhere in the system. This is the case if f{^*) < 0. Therefore, the region of 
parameters being considered is defined by 



k > 



+sinh2Ci 



The exponential WKB approximation in this case reads 

bi 



VWx 



: exp 



VWx 



: exp 



(42) 



(43) 



with 



H^x) = V^li^), (44) 
and hi and 62 arbitrary constants. Imposing the boundary conditions (pO) and (p^) leads to the relationship 



MO)^tanh /^^de.Mex) = ^^^J^ 
Jo cosh" £_* 

For a thermal wall at .t = 0, the above equation is substituted by 

cosh / d^x h{£,x) = 0, 
Jo 



(45) 



(46) 



and, therefore, there is no WKB solution in this range of parameters. 

c) The function f{£,x) changes sign in the interval < < It must be /(O) < and /(^^) > 0. This corresponds 
to the k interval 



+ sinh 2^* ^ ^ ^ 2^* + sinh 2C 



2^2^* cosh^ 2^* 



2^/2e^ 



(47) 



Since /(^k) exhibits in this case a zero in the integration range, we have to consider separately the regions < o- 
and ^x > CL, where a is the turning point, i.e. f{a) — 0. In the former region, the WKB solution is given by an 
expression of the form (^), while in the latter it has the form (^). A global solution is constructed by matching 
both WKB approximations through the connection formulas, expressing the connection between the oscillatory and 
exponential behaviors. Since it is easily seen that the turning point is a simple (first-order) zero, a standard application 
of the theory suffices ||l5| . Then, imposing the boundary conditions gives the following equation to be verified by the 
solutions in this range: 
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TT \ 3 tanh - (0) cosh^ Q 1 



3tanhC* + /i3(0)cosh2^* 



■ exp 



-2 / 

'0 



where 

J a 

The resuh for a thermal wall at x = is given in this case by 



tan ~ ^) = ^ "^XP -2 / hi^^ 

L ■'0 



(48) 



(49) 



(50) 



Once we have derived the equations determining the solutions to the boundary value problem, the strategy to 
build up the marginal stability curve is as follows. Given a value of G, the first question is to see whether there 
is a bifurcation from the one-dimensional solution, i.e. whether Eq. (pO) has a solution. In case the answer to this 
question is affirmative, the smallest value of f * for which there is solution gives the point of the marginal stability 
curve corresponding to that value of For larger values of the one-dimensional solution is not stable. Suppose 
a solution of Eq. (|3^) exists for a given wavenumber k. This value is compatible with many (infinite) values of 
the smallest one corresponding to the choice g = 1 in the relationship between the possible values of k and ^* (see 
below Eq. (^)). Moreover, the larger the value of k, the smaller the value of ^* associated. The conclusion is that 
the stability curve is determined by the largest value of k for which the problem has solution. From a practical 
point of view, we have to start by looking for solutions belonging to case b discussed above. The values ^* = n/k 
obtained from the solutions of Eq. (^5|) are the critical values, i.e., they define points on the marginal stability curve. 
Nevertheless, there is no solution for every value of ^* belonging to the fc-interval defining case b. More precisely, 
there is no solution at all for a thermal wall, while it is a simple matter to show that the existence of solution in this 
region for a vibrating wall requires that 



tanh 



< 



3 tanh 
cosh^ ^* 



where 



/i(i)(C.) = lim /i(C.). 

k—>-ki 



(51) 



(52) 



Equation (|l|) gives f * ~ 0.555. Therefore, solutions for larger values of ^* must be found, if they exist, in a different 
region of values of fc, namely in that considered in case c. With this procedure, it is an easy task to find the largest 
value of k (smallest value of ^*) for which the eigenvalue problem defined by Eq. (^0|) has solution for each value of 
f*, in the WKB approximation. The marginal stability curve for a vibrating wall obtained in this way is given in 
Fig. ^ where we have plotted the critical values Ac of the asymmetry parameter A = Ly/L^ = ^*/^*, as a function 
of Q . The solid line is the WKB solution corresponding to the so-called case 6, while the dashed line is from case c. 
Note that the solutions from both regions of parameters match rather smoothly at ^* ~ 0.555. Moreover, the critical 
asymmetry Ac is a monotonic decreasing function of growing very fast when Q tends to zero, as it was expected 
since the transversal inhomogeneities must disappear in the elastic limit. 

For a thermal wall, we have obtained that there is no solution in region b. Therefore, we must search for solutions 
belonging to the range c of parameters. Analysis of Eq. (^ ) leads to the result that the equation has no solution for 
^* ^ 1.181. Moreover, Eq. ( |4l| ) has no solution in that region either, implying that, in the WKB approximation, the 
transition to the transversally inhomogeneous state requires, in the case of system driven by a thermal wall, that the 
inelasticity and the number of monolayers at rest be not small. The dotted line in Fig. |^ shows the marginal stability 
curve for a thermal wall. In the limit of large ^* , the curve overlaps with the one for a vibrated system. 



IV. SIMULATION RESULTS 



To test the above theoretical results and also to investigate the nature of the predicted symmetry breaking, two 
more fundamental descriptions of the system, via the direct simulation Monte Carlo (DSMC) method and molecular 
dynamics (MD) simulation, respectively, have been considered. Although both are based on a dynamical simulation of 
the system, their nature is rather different. The DSMC method provides an algorithm to obtain numerical solutions 
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of the Boltzmann equation for given initial and boundary conditions Therefore, it rehes on the vahdity of a 

kinetic theory description, namely the one given by the Boltzmann equation. On the other hand, no hydrodynamic 
approximations are introduced into the description, so that the validity of a hydrodynamic level of description, as 
provided for instance by the Navier-Stokes equations, is not taken for granted in this approach. 

The MD simulations follow the motion of the particles of the system as a sequence of free motions and binary 
collisions [Q, i.e. by direct application of Newton's equations of motion. Therefore, MD provides the more basic 
description of the evolution of the system. In this context, it is important to stress that although the DSMC method 
also uses "particles" at a formal level, they are not real particles, but fictitious ones, which are employed to mimic the 
ideal dynamics described by the Boltzmann equation. It is precisely the ideal nature of the particles in the DSMC 
method what allows a very high numerical accuracy, since the number of particles used do not affect at all the physical 
state being simulated. In particular, the above number is not related with the actual density of the system. 

In the following, we will report simulation results obtained for the system studied in the previous sections, restricting 
ourselves to the vibrating wall. We have used the two methods, DSMC and MD, because they complement one 
another. For instance, comparison of the results obtained by both methods provides a test of the validity of the 
inelastic Boltzmann equation to describe a system under the physical conditions used in the MD simulations. Since 



the one-dimensional state has been discussed in detail elsewhere [10;!^, attention will be focused here on the behavior 
of the system in the vicinity of the bifurcation. 

In the MD simulations, the density plays of course a relevant role. As the interest here is in the low density limit, 
small values of the surface fraction v = Ni^a^ / AL^Ly, typically of the order of 10~^, have been used. For fixed 
given values of a and L^, a set of simulations have been run corresponding to different widths Ly of the system, 
starting from a small enough value. This means that each set of simulations corresponds to the same value of as 
defined in Eq. (p8|), while they differ in the asymmetry parameter A. The simulations of the Boltzmann equation 
have been carried out with a similar systematic, the main difference being that the density plays no role in them. 
The other parameter needed to specify the simulations is the velocity of the vibrating wall Vb- We have verified that, 
in agreement with the theory developed in the previous sections, the formation of transversal inhomogeneities is not 
altered by modifying this velocity, as long as it is large enough to fluidize the complete system. 

All the simulations started from a spatially homogeneous configuration, with a Gaussian velocity distribution, 
corresponding to an arbitrary temperature. The simulation is then followed until the system reaches a steady state, 
in which the several monitored properties of the system (mean kinetic energy, density fluctuations, and hydrodynamic 
profiles) become time independent. Once the system is in the steady state, all statistical averages of interest are 
accumulated. For the purpose here, the density and temperature profiles provide the relevant information. Let us 
describe what is observed in a set of simulations, namely we are going to present DSMC data from systems with 
a = 0.95 and — lOA. This corresponds to ^* = 1.015. In Fig. |^ the steady two-dimensional density profile for 
Ly = 20A (A = 2) is shown in a three-dimensional plot. It is seen that no gradients in the y-direction are present, 
i.e. the system is in the one-dimensional reference state. The transversal homogeneity appears even clearer in Fig. |^, 
where the density is plotted as a function of y for several fixed values of x. 

When the width Ly is increased keeping all the other parameters fixed, a critical value shows up such that gradients 
in the y-direction spontaneously develop in the system for larger widths. An example of this is given in Fig. ^ where 
the density surface for the same parameters as in Fig. |^, except that now Ly — 26. 5A, is plotted. A gradient in the 
direction parallel to the vibrating wall is clearly identified, becoming more pronounced next to the opposite wall, 
as illustrated in Fig. The density gradients in this case are relatively small, being the maximum variation of the 
density in the y direction of the order of 5%. In fact, in all the simulations, both by DSMC and MD, it has been found 
that when increasing Lj,, a continuous transition from the one-dimensional state to a state with weak inhomogeneities 
in the transversal direction occurs. Moreover, the transversal density profile near the transition exhibits, as in the 
case of Fig. ||, a wavelength equal to twice the width of the system. In the language used in Sec. Ill, what is observed 
is the perturbation with q = 1 (k = tt/^*), consistently with our theoretical analysis. 

If the width of the system is increased further on, the transversal inhomogeneities grow very fast and, of course, 
the results from the linear marginal stability analysis do not apply. Simulations show that the density becomes very 
large in one of the corners of the system, away from the vibrating wall. The rapid increase of the gradients and the 
large value of the density in a localized region of the system, lead to conclude that in this regime a cluster of particles 
IS formed §. Of course, for such a region of parameters the hydrodynamic profiles obtained from DSMC and MD are 
quite different, as the former considers the particles as points while the latter assigns them a finite diameter, implying 
that the density is bounded by the close-packing value. 

Once the appearance of transversal inhomogeneities has been observed by visual inspection, it is desirable to have 
a "quantitative" criterium to establish whether the system is or is not transversally homogeneous. This is equivalent 
to identifying an order parameter to characterize the transition. Since there may be gradients in both directions, the 
identification of such a parameter is not at all a trivial task. The one we have chosen is defined as follows. First, we 
introduce the dimensionless quantity Pxiu) by 
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If the system is transversally homogeneous, Px{y) is independent of both x and y, and equal to zero. When transversal 
inhomogeneities are present, it depends of course on y but, as it can be guessed from Figs. ^ and |^, it also depends 
on X. Nevertheless, the simulation data indicate that this dependence is rather weak, at least near the transition. 
As an example, in Fig. |^ the function Pxiy) has been plotted for the same system as in Figs. ^ and ||. The different 
lines correspond to four different values of x, equally separated, in the interval [3X3^/4, L^]- This is the region where 
the transversal gradients are larger. From the figure it follows that the x dependence is practically scaled out in 
the definition of Px{y)- It must be mentioned, however, that if the whole range of variation of x is considered, some 
dependence on x shows up. In any case, the departure from zero of the average value of p-xiy), piy), over a certain x 
interval, next to the vibrating wall and not too large to avoid x dependence, provides a good criterium to distinguish 
transversally inhomogeneous systems from homogeneous ones. The results to be discussed in the following have been 
obtained using the interval [3Lx/4:, L^]- 

The above discussion, the theoretical analysis, and the numerical results, as those illustrated in Fig. |^, suggest that 
a good order parameter may be given by the absolute value of the first Fourier component, of piy). In fact, 
we have verified that the absolute value of the Fourier transformed of this quantity exhibits an abrupt maximum for 
the first component when transversal gradients begin to build up in the system. The behavior of |/i| as a function 
of the asymmetry A in the vicinity of the transition point, is shown in Fig. Q for the same parameters considered 
in Figs. I^joj First of all, it must be noted that |/i| varies in a continuous way through the transition, although it 
grows very fast when one goes into the inhomogeneous region. This is the typical behavior of the order parameter of 
a nonequilibrium second order phase transition p^ . The continuous character of |/i| introduces some arbitrariness in 
the determination of the transition point Ac. We have made the choice, somewhat arbitrarily but consistently, that 
the transition takes place when |/i| becomes an order of magnitude larger than its typical value in the reference state, 
which is determined by the noise level. For the systems used in DSMC, this latter value is of the order of 10~^, so 
that a system has been considered as transversally inhomogeneous when |/i| ~ 10~^, which implies deviations from 
homogeneity of the order of 1%. This leads in the case of Fig. |^ to an estimation Ac = 2.5 ± 0.1 for the critical 
asymmetry. 

By changing the initial parameters of the system and repeating the above procedure, Ac has been computed for 
different values of The results from DSMC are represented by the filled symbols in Fig. |l} The agreement between 
the theoretical predictions and the simulation data is rather good, although a systematic deviation appears, larger the 
smaller When valuing the comparison, it must be taken into account that the simulations only provide an upper 
bound for Ac. When the system is very close to the transition point, the time required to go from the transversally 
homogeneous state to the inhomogeneous one may be too large to observe the transition during the simulation time. 
In any case, it is fair to say that the hydrodynamic equations and the WKB approximation provide an accurate 
description of what is observed in the simulations. 

The results discussed up to this point were obtained from DSMC. Just to illustrate how MD simulations lead to 
an equivalent scenario, we present next some results for a set of MD simulations with a = 0.925, v = 10~^, and 
Lx = lOOfj. For these values, it is ^* = 0.281. In Fig. || a three-dimensional plot of the density profile is shown 
for A = 6.4. The system exhibits gradients in the transversal direction, increasing again as we move away from 
the vibrating wall. However, the density gradients are very small, being the maximum deviation from homogeneity 
of the order of 10%. That means that the system is close to the transition, probably in the region where a linear 
approximation around the steady state still provides an accurate description. 

When the asymmetry is increased further on, gradients in the perpendicular direction become sharper, and a state 
with a very peaked density is observed. This is illustrated in Fig. ^ for A = 7. In Fig. the quantity p is shown for 
different simulations belonging to the set we are considering, i.e. they only differ in the value of A. The continuous 
transition from the reference state to the transversally asymmetric one is clearly observed, as well as the dramatic 
increase of the transversal gradients when the system gets well inside the unstable region. In conclusion, MD results 
are in full qualitative agreement with the DSMC ones. Even more, the critical values of the asymmetry parameters 
obtained from MD are in very good quantitative agreement with those following from DSMC, as seen in Fig. |l|, where 
they have been represented by the open symbols. This provides a proof of the validity of the kinetic theory description 
as given by the Boltzmann equation for dilute inelastic gases. 

It is worth to emphasize that points in Fig. |l| were obtained from DSMC and MD simulations by changing the 
values of a and aN/Ly to sample different values of The fact that Ac obtained in this way varies smoothly with 
^* supports the theoretical prediction that the dependence on the different parameters of Ac occurs through it. This 
has also been ratified by considering two sets of simulations having different values of both a and aN/Ly, but leading 
to the same value of The same critical asymmetry parameter was found, supporting the scaling predicted by the 
theory. 



V. DISCUSSION AND CONCLUDING REMARKS 



In this paper, the spontaneous transversal symmetry breaking in a vibrated granular fluid predicted by Livne et al. 
has been further investigated for low density systems. It has been shown that the transition takes places both in 
MD simulations and in systems described by the Boltzmann equation. Moreover, there is a reasonable good agreement 
between the theoretical predictions, following from a marginal stability analysis of the hydrodynamic description of 
the system, and the results from MD and also from the numerical solution of the Boltzmann equation by the DSMC 
method. This refers to the values of the critical asymmetry as a function of the control parameter, and also to the 
shape of the hydrodynamic profiles in the vicinity of the symmetry breaking. 

To characterize the transition, an order parameter quantifying the initial set up of transversal inhomogeneities has 
been introduced. In terms of this parameter, the transition presents the features of a second order nonequilibrium 
phase transition. In this context, it is worth mentioning that no sub-critical bifurcations have been observed in the 
simulations. Moreover, for states well inside the instability curve, the density profile exhibits a characteristic nonlinear 
A/2 shape. On the other hand, in ref. 0, nonlinear two-dimensional states inside the lineal stability region were found 
at high densities from the numerical solutions of the hydrodynamic equation. Of course, there is not any contradiction 
in this, since our analysis has been restricted to low density gases. In ref. |^ an analytical expression for the marginal 
stability curve in a certain limit is given. In our notation, the limit considered is ^* 1, and the expression reads 
A ~ 1.6/^*^, where we have neglected subleading terms in the density. An asymptotic analysis of the WKB results 
in this paper, namely of Eq. (45) leads to Ac ~ 1.63/^*, i.e. a qualitatively different behavior. Given the asymptotic 
character of the region where these expressions are derived, it is hard to discriminate between the two results from 
the simulation data. 

The work reported here shows once again the generality of spontaneous symmetry breaking phenomena in granular 
fiuids. They occur in isolated as well as in driven granular systems, in dense and dilute flows, with and without 
gravitational field acting on the particles, i.e., they appear as quite a ubiquitous effect. Interestingly, all the indications 
point out that they always have a collective origin, which is fully captured by a hydrodynamic description. 

Although it may be thought that the spontaneous symmetry breaking discussed in this paper must be closely related 
with the clustering instability exhibited by a freely evolving granular gas |^ , we believe this relationship deserves more 
work. The information we have up to now about each of the two phenomena is not the same. We know the mechanism 
responsible for the clustering instability (the growth of the shear mode relative to the granular temperature) , but not 
the final state attracting the system. On the other hand, the existence of a transversally inhomogeneous steady state 
has been established for vibrated systems, but the detailed mechanism responsible for the development of instabilities 
from the reference state is not known. 
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FIG. 1. Marginal stability curve Ac(^^). The solid and dashed lines are the WKB predictions for a system driven by a 
vibrating wall, while the dotted line is for a system with an isothermal wall. Also shown are the results from DSMC (filled 
symbols) and MD (open symbols) simulations. For a given value of Q, the system exhibits transversal inhomogeneities in the 
steady state above the marginal curve. 

FIG. 2. Three-dimensional plot of the stationary density profile obtained by the DSMC method for a system with A = 2 
and Q — 1.015. The density is normalized by its average value n, and the lengths with the average mean free path A. The 
system does not exhibit gradients in the j/-direction, i.e. it is in the reference state (bellow the marginal stability curve). 

FIG. 3. Density profiles as a function of y for several fixed values of x, for the same system as in Fig. |l] The curves 
correspond, from bottom to top, to a; = Lx/4:, Lx/2, 3Lx/4:, and L^- 

FIG. 4. Three-dimensional plot of the density profile obtained by the DSMC method for the same system as in Fig. ^ with 
the only difference that now it is hy = 26. 5A. The system is now above the marginal stability curve and transversal gradients 
are clearly observed. 

FIG. 5. Density profiles as a function of y for several fixed values of x for the same system as in Fig. ^. The curves correspond, 
from bottom to top, to a; = Lx/4:, La:/2, SL^/i , and Lx- 

FIG. 6. DSMC results for the dimensionless function Px{y), defined by Eq. (p^, for the same system as in Figs. |^ and ^. 
The different curves correspond to equally spaced values of x in the interval [3La;/4, L^]- 



FIG. 7. DSMC results for dimensionless order parameter [/ij as a function of the asymmetry A for a system with Q ~ 1-015 
(the dashed line has been included as a guide for the eye). A second order nonequilibrium bifurcation is clearly identified. 

FIG. 8. Three-dimensional plot of the stationary density profile obtained by MD simulation. The particles are disks of 
diameter a, the area fraction is i/ — 10~^, the coefficient of restitution a = 0.925, = lOOcr, and A = 6.4. Small transversal 
gradients are present, and the system is in the vicinity of the bifurcation. 

FIG. 9. The same as in Fig. ^ but with A = 7. Quite large transversal gradients growing in the x direction are identified. 
The system is above the marginal stability curve. 

FIG. 10. MD results for the dimensionless function p{y), for several values of the asymmetry A, as indicated in the figure. 
All the other parameters of the system are the same as in Figs. ^ and |^. The transition to a state inhomogeneous in the y 
direction is clearly identified. 
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